/*
 * AqueousActivity.cpp
 *
 *  Created on: 19 Aug 2011
 *      Author: Allan
 */

#include "AqueousActivity.h"

// C++ includes
using namespace std::placeholders;

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include <Geochemistry/Activity/AqueousModels/HkfModel.h>
#include <Geochemistry/Activity/AqueousModels/SetschenowModel.h>
#include <Geochemistry/Activity/AqueousModels/DrummondCO2Model.h>
#include <Geochemistry/Activity/AqueousModels/RumpfCO2Model.h>
#include <Geochemistry/ElectrolyteSolution.h>
#include <Numerics/LinearInterpolator.h>
#include <Utils/Algorithms.h>

AqueousActivity::AqueousActivity()
{}

AqueousActivity::AqueousActivity(const string& species, const vector<string>& phase_species, const ActivityCoefficientFunction& activity_coefficient) :
activity_coefficient(activity_coefficient)
{
	index_species = IndexOf(species, phase_species);
	
	index_water = IndexOf("H2O", phase_species);
}

AqueousActivity::~AqueousActivity()
{}

const double AqueousActivity::activity(double T, double P, const VectorXd& n) const
{
	const double molality_species = 55.508 * n[index_species]/n[index_water];
	
	return molality_species * activity_coefficient(T, P, n);
}

const double AqueousActivity::activityCoefficient(double T, double P, const VectorXd& n) const
{
	return activity_coefficient(T, P, n);
}

const VectorXd AqueousActivity::activityDerivatives(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	VectorXd derivatives = VectorXd::Zero(n.rows());
	
	derivatives[index_species] = a[index_species]/n[index_species];
	
	return derivatives;
}

const AqueousActivity AqueousActivityIdealModel(const string& species, const ElectrolyteSolution& electrolyte_solution)
{
	const Index index_water = electrolyte_solution.GetWaterIndex();
	
	auto water_molar_fraction_function = [=](double T, double P, const VectorXd& n) 
	{
		return n[index_water] / n.sum(); 
	};
	
	AqueousActivity aqueous_activity(species, electrolyte_solution.GetSpecies(), water_molar_fraction_function);
	
	return aqueous_activity;
}

const AqueousActivity AqueousActivityHKFIonicModel(const string& ion, const ElectrolyteSolution& electrolyte_solution)
{
	const vector<string> ions = electrolyte_solution.GetIons();
	
	vector<double> effective_electrostatic_radii_ions;
	
	BOOST_FOREACH(const string& iter, ions)
		effective_electrostatic_radii_ions.push_back(EffectiveElectrostaticRadius(iter));
	
	const Index local_index_ion = IndexOf(ion, ions);
	
	ActivityCoefficientFunction activity_coefficient_function = 
		std::bind(HkfIonicActivityCoefficient, electrolyte_solution, local_index_ion, effective_electrostatic_radii_ions, _1, _2, _3);
	
	AqueousActivity aqueous_activity(ion, electrolyte_solution.GetSpecies(), activity_coefficient_function);
	
	return aqueous_activity;
}

const AqueousActivity AqueousActivityHKFWaterModel(const ElectrolyteSolution& electrolyte_solution)
{
	const vector<string> ions = electrolyte_solution.GetIons();
	
	vector<double> effective_electrostatic_radii_ions;
	
	BOOST_FOREACH(const string& iter, ions)
		effective_electrostatic_radii_ions.push_back(EffectiveElectrostaticRadius(iter));
	
	ActivityCoefficientFunction activity_coefficient_function = 
		std::bind(HkfWaterActivityCoefficient, electrolyte_solution, effective_electrostatic_radii_ions, _1, _2, _3);
	
	AqueousActivity aqueous_activity("H2O", electrolyte_solution.GetSpecies(), activity_coefficient_function);
	
	return aqueous_activity;
}

const AqueousActivity AqueousActivitySetschenowModel(const string& species, const ElectrolyteSolution& electrolyte_solution)
{
	return AqueousActivitySetschenowModel(species, electrolyte_solution, 0.0);
}

const AqueousActivity AqueousActivitySetschenowModel(const string& species, const ElectrolyteSolution& electrolyte_solution, double b)
{
	return AqueousActivitySetschenowModel(species, electrolyte_solution, {b}, {25.0}); // Assume the value b is given at 25 °C
}

const AqueousActivity AqueousActivitySetschenowModel(const string& species, const ElectrolyteSolution& electrolyte_solution, const vector<double>& bi, const vector<double>& Ti)
{
	LinearInterpolator b(Ti, bi);
	
	ActivityCoefficientFunction activity_coefficient_function = 
		std::bind(SetschenowActivityCoefficient, electrolyte_solution, _1, _2, _3, b);
	
	AqueousActivity aqueous_activity(species, electrolyte_solution.GetSpecies(), activity_coefficient_function);
	
	return aqueous_activity;
}

const AqueousActivity AqueousActivityDrummondCO2Model(const ElectrolyteSolution& electrolyte_solution)
{
	ActivityCoefficientFunction activity_coefficient_function = 
		std::bind(DrummondCO2ActivityCoefficient, electrolyte_solution, _1, _2, _3);
	
	AqueousActivity aqueous_activity("CO2(aq)", electrolyte_solution.GetSpecies(), activity_coefficient_function);
	
	return aqueous_activity;
}

const AqueousActivity AqueousActivityRumpfCO2Model(const ElectrolyteSolution& electrolyte_solution)
{
	ActivityCoefficientFunction activity_coefficient_function = 
		std::bind(RumpfCO2ActivityCoefficient, electrolyte_solution, _1, _2, _3);
	
	AqueousActivity aqueous_activity("CO2(aq)", electrolyte_solution.GetSpecies(), activity_coefficient_function);
	
	return aqueous_activity;
}
